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ABSTRACT 

The in-plane edge impact of composite plates, with or without a pro- 
tection strip, is investigated in this work. A computational analysis 
based on the Fast Fourier Transform technique is presented. The particular 
application of the present method is in the understanding of the foreign 
object damage problem of composite fan blades. However, the method is 
completely general, and may be applied to the study of other stress wave 
propagation problems in a half space. Results indicate that for the pro- 
tective strip to be effective in reducing impact stresses in the composite 
the thickness must be equal or greater than the impact contact dimension. 
Also large interface shear stresses at the strip - composite boundary can 
be induced under impact. 
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INTRODUCTION 

This report is part of a continuing effort by NASA to understand the basic 
mechanics of foreign object impact of composite materials of particular interest 
are damage resistant designs of jet engine fan blades under hail or bird 

impact. In previous reports the central or normal impact response of composite 

plates was examined [1], In this report the mechanics of edge impact of 

composite plates are examined. This is schematically illustrated in Figure 1. 
The basic approach to the study of impact of composite plates in this 

program has been to examine the stress waves generated by the impact forces. 

For central impact of plates it has been shown that in addition to wave pro- 
pagation across the plate thickness, bending and extensional waves propagate 
away from the impact site. The stresses associated with these waves have been 
studied without considering the effect of boundaries such as the free or clamp- 

led edges of a fan blade. This simplification has been made on the premise that 

-4 

for short impact times e.g. less than 10 sec. few edge reflections have 
taken place, and that the highest stresses occur at the impact site. However, 
for edge impact, the boundary conditions greatly affect the nature of the wave 
mechanics. 

Edge waves in solids have been studied extensively in seismology. The 
principal phenomenon is the entrapment of wave energy in a layer near the sur- 
face. This surface wave is known as a Rayleigh wave and travels at a velocity 
below the shear velocity for isotropic solids. For plates, however, two types 
of edge waves can occur as shown in Figure 2a. For impact transverse to the 
plate, flexural edge waves can occur. For in-plane, Rayleigh type 
edge waves are generated. In this report only in-plane edge waves will be 
discussed. 
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Wave type solutions to the equations of elastodynaroics of an orthotropic 

plate which exponentially decay away from the edge (X direction) and pro- 

3 

pagate along the edge (X direction) can be found provided the edge wave 

1 

velocity satisfies the equation (Reference 9) . 


pv 2 + 


C - pv 2 
55 

C C (C - 
33 55 11 


pv 2 ) 


1/2 2 

[C - C (C - pv 2 )] = 0 
13 33 11 


where P is the massdensity of the plate and are the effective plate 
elastic constants (denoted by / C V in Reference 1). It can be shown that 

ij 

one real root lies in the. interval 

0 < pv 2 < C 

55 

Thus, the edge or Rayleigh wave speed is less than the shear speed in this 
1/2 

direction [C /p ] 

55 

Changing the layup angle will affect the elastic constants C and hence, 

ij 

change the value of the edge wave velocity. The results of this calculation 

are shown in Figure where the C are obtained from Reference 7 . The 

ij o 

edge wave speed seems to obtain a maximum between +_ 15 and +_ 30° layup angle which 
is below the extentional wave speeds labelled "dilational" and ’'shear 11 
in Figure 2b and which is greater than the bending wave velocity. 
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In order to prevent damage to composite fan blades under foreign object 
impact leading edge protection has been used. (See e.g. Ref. (8)). This usually 
consists of a strip of metal attached to the leading edge of the fan blade. To 
model the effects of this impact protection strip, the in-plane edge impact 
of an anisotropic plate, with a beam- strip attached to the impact edge, has 
been studied (see Figure 1). It will be shown later that the strip will decrease the 
tensile stress along the edge while producing shear stress between the strip 
and the plate edge. 
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I. BASIC EQUATIONS FOR EDGE LOADING ON ANISOTROPIC PLATE 

Let the plate be the half-space x 3 > 0. The equations of motion 
are given by [1] as: 


C u + C u + (C +C)u = p Uj 

11 1,11 55 1,33 13 55 3,13 ,tt 


C U + C U + (C +C)u 

33 3,33 55 3,11 13 55 1,13 


P u 


3,tt 


( 1 . 1 ) 


here we have employed C , C , C ... to denote C , C , C ... of 

11 13 33 11 13 33 

[1], The in-plane motion is assumed to be independent of the bending defor- 
mation. With the loading condition shown in Figure 1, the boundary conditions 
are (without protection strip) : 


t (x ,0,t) = C (u 

13 1 5 5 1 


, 3 


+ 


U 

3,1 


) 


x =0 

3 


0 


( 1 . 2 ) 


t (x ,0,t) 
3 3 1 


[C u + C 

3 3 3 , 3 1 3 


u .1 I = p f (X ) g(t) 

1,1 ] X =0 1 

3 


The following nondimens ional parameters are used; 
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C * = C /C , C* = C /C , C* = C /C , C* = C /C , P* = p/C 

1 ] 11 66 13 13 66 33 33 66 55 55 66 0 66 

(1.3) 


u* = u JSL, u* = u /l, x* = x /£, x* = x /£, t* = t /c~~ /p 'fi (1.4) 

11 33 11 33 66 


i is a length parameter. C is a typical elastic constant. In what 

6 6 

follows the equations will be assumed to be nondimens ionali zed using (1.3) 


To obtain the solution, we employ transform methods. Define: 


F(f) = 


e“ lk i x i f(xi) dxi 


as the Fourier transform of f(x) where we assume that 


f(x ) = j> f(x ) = 0 

1 ^ x. 1 

_oo 1 _oo 


F(3 f) = ik F (f ) , 

*1 


F(3 2 f) = -k 2 F (f ) 


as the Laplace transform. The initial conditions are assumed to be. 
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u (x ,x ,0) = 3 u (x,,x ,0) = 0 
113 a* i 1 3 


(1.5) 


u (x .x ,0) = 3 u (x ,x ,0) = 0 

3 13 Yt 3 1 3 


Letting £jP(u j] = u (k,x ,s),f[F(u )}= u (k,x 
1 1 3 3 3 3 

transformed equation: 


,s), we have the following 


C (-k 2 ) u + C u + (C +C ) (ik ) u = +s 2 u 

11 1 1 551,33 1355 1 3,3 1 


( 1 . 6 ) 


C u + C (-k 2 ) u + (C +C ) (ik ) u = +s z u 

333,33 55 1 3 1355 1 1,3 3 


2 " 


and at x = 0 , the transformed boundary conditions become 

3 


u + ik u = 0 

1,3 13 

(1.7) 


C u + C ik u = -p f(k )g(s) 

333,3 13 11 ° 1 


Since we are expecting surface wave propagation, we seek the solution of 


(1.6), (1.7) in the forms: 


u 

1 

= 

'O 

e -p(k,s)x 3 


U 

3 


* 

2 


(real (p) > 0) 


0 . 8 ) 



7 


Therefore, the equations for 4> , $ are: 

l 2 


s 2 - C k 2 + p 2 C -ik p(C +C ) 

11 1 r 55 I 13 55 


ik p(C +C ) 

1 13 5 5 


- s 2 - k 2 C + p 2 C 

15 5 


' 

-e-l 


l 


4> 


2 


= 0 


( 1 . 9 ) 


or for a non-trivial solution, 


det = C C p* ♦ [C (~s 2 -C k 2 ) + C (-s 2 -C k 2 ) + k z (C +C ) z ] p z 

3355 33 111 55 551 11355 


+ (s 2 +C k 2 ) (s 2 +C k 2 ) = 0 

111 551 


(1.10) 


We will choose the p’s with positive real parts to ensure the decay in x 

3 

direction of the surface wave. Let the solutions be p = p ,p , therefore, 
we have 


Cl) 

P * E C^k^s) , <P 2 = 


i [-s -C k 2 ♦ C p 2 ] 

I 1 1 5 5 1 

k p (C + c ) 

II 13 55 


<f>^ £ Ip c 

1 3 11 


(l.H) 


2 ^) 

p = p : ^ ^ C (k ,s), <p = C = 

*21 212 32.2 


i [ - s 2 - C k 2 + C p z ] 

_ 1 1 1 5 5 2 

k p (C +C ) 

12 13 55 


( 2 ) 


Therefore, the displacements have the forms: 


U 

1 

u 

II 

i 

e ^ i X 3 + C ^(k ^,s) 

1 

U 

L 


* 

_ 3 il 


V 

^ 3 2_ 


( 1 . 12 ) 
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C , C are determined from the boundary conditions (1.7), or 
1 2 


- p . 

+ il 

< *P , 

1 31 

- p 

- p . 

C 

3 3 

ijj + ik C 

3 1 113 

- p 


1 32 


-P 0 H 


(1.13) 


Here, the determinant, A(p ,s), given by, 

1 


A - (ik \p - p ) (ik C 13 - p C ip ) 

1 31 1 1 2 33 32 


- (ik ip - p ) (ik 

1 32 r 2 1 


c 

1 3 


- p C ^ ) 

*1 33 31 


(1.14) 


must be non-zero to ensure a solution. (A=0 gives the Rayleigh poles of 
the system, which correspond to a free surface). Therefore, 


C 

i 



4> 

3 2 


- p 2 ] fg p 0 


c 

2 



ik 




f g P, 


(1.15) 


and the physical displacements 
inverting their transforms. 


u (x ,x ,t), u (x , x ,t) 
113 3 13 


are obtained by 
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From the stress-strain relations 


t 

1 3 


C 

55 


(U . +U 
1,3 3 


), 


t = C u + C u , 

11 111,1 13 3,3 


t=Cu+Cu (1.16) 

33 33 3,3 13 1,1 


we obtain the transforms of the stresses 


t 

1 3 



U-P 


+ ik 

l 


) C 

3 1 1 


-p X 

e r i 3 


+ (-p +ik ip ) C e ^ 2 X 3 } 

2 13 2 2 


t = (C ik -C p ip ) C e' P i X 3 (1.17) 

11 11 1 .13 13 1 1 


(C ik -p ip C ) C e' P 2 X 3 

11 1 2 32 13 2 


t 

3 3 


(-p lp C +ik C ) 

13133 113 


-p X 
e r i 3 


+ (-p ip C +C ik ) C e P 2 X 3 , 

2 3 2 3 3 1 3 1 2 
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and the physical stresses can be obtained by inversion. 
The particular forcing function employed is 


2 

f (V e(t) = [1 - (li) ] sin II 

l o 


0 < t < x J x I < a 

0 i — 



( 1 . 18 ) 


where x* = l/J~C 7p~^ 

0 66 


Here a is a length measuring the impact area, and t q is the contact 
time [1]. The transform of this particular forcing function, in non-dimen- 
sional form becomes, 


f g = 


kj (a/ A) 


ft 


cosk^C^) f sin k (|) 


VT*Cl + e _ST S) 


+ T* 2 s 2 


( 1 . 19 ) 
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II. EDGE IMPACT OF PLATE WITH EDGE PROTECTION 

To prevent failure of composite fan blades under impact forces, leading 
edge protective strips have been employed. In practice, these strips of 
stainless steel are wrapped around the leading edge. To model this device, 
we consider a beam bonded to the edge of an anisotropic plate (Figure l)* 

The effect of the beam will be to thwart the force of impact, thereby 
decreasing the normal stresses in the composite. However, we shall show 
that with such a reduction in normal stress , sizeable interface shear stresses 
can be induced. 

With the introduction of a beam of thickness b on the edge of the composite 
plate, the Rayleigh wave behavior will depend on the ratio of the wavelength 
to thickness ratio of each Fourier component in the direction. Thus one 
should expect the Rayleigh wave speed to vary with b/a, the thickness to 
impact footprint ratio. In addition the Rayleigh wave will become distorted 
as it propagates. 

To solve the edge strip problem the solution in the composite plate 
follows the same procedure as the no- strip case except for the boundary 
conditions on the edge. In place of the zero stress conditions on the edge 
we relate the edge stresses t^» t^ to the motion of the beam strip. If 
one considers a small element of the beam-strip along the x^ direction, the 
momentum balance equations in the x^, x^ directions become, (for a plate 
of unit thickness) 

pb iiLL = Eb iiE + b 13 

at 2 


( 2 . 1 ) 
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pb 3 2 W = - El 3 4 W + I b 3 4 W + b at,- + t__ + pfg (2.2) 

p ~2 w o 

3t 2 3x 4 ! 3x 2 3t 2 3X]_ 

4 

In these equations U, W are the x^, displacements of the beam element at the 
half thickness, and t^, t^ are "the interface stresses. 

We choose the compatibility conditions between the beam and plate dis- 
placements 

W = u^, on x^ = 0. (2.3) 


U 




on x^ = 0 


(2.4) 


In the above equations b is the depth of the strip, E, I, Ip are 
respectively the Young's modulus, moment of inertia and rotary inertia. Also 
p f(t) g(Xj) is the edge loading now applied to the outer protective strip 
surface. 

The equations for the plate remain as in the free edge case and a 
solution is obtained by taking a Laplace transform on time and a Fourier 
transform on the space variable x^ . With nondimensionalization the solution 
in the plate is assumed in the form of (1.12). The transform of the plate 
displacements are 

i 

e" P 2 X 3 (2.5) 

where p^, are defined in (1.10) and are given in (1.11). C^, 

are determined from the edge boundary conditions. However, in place of the 
free edge conditions (1.2) we use the equations of motion for the strip (2.1), 


u. 


u, 


= C. 


e" P l X 3 + C, 


31 


32 



13 - 


(2.2). C^, C 2 are then solutions to the algebraic equations 


G. 

H, 


— — 
C, 


0 

1 

1 


1 


G 

H 


C 2 



- 



L — 


— ^ 


where 


( 2 . 6 ) 


G = -k 2 Ew + p C - pws 2 - 4; (ik 2 w 2 E/2 + ikC + pw 2 s z ik/2) 

l l 55 31 55 


H = -k 2 Ew + p C - pws 2 - ip (ik 3 w 2 E/2 + ikC + pw 2 s 2 ik/2) 
l r 2 55 32 55 


G = "ikC + ikwp C/2 + ( p C + Ew 2 k 4 /12 + wk 2 C /2 

13 55 3 1 M 3 3 55 


+ k 2 s 2 pw 3 /|2 + pws 2 ) 


H 


= -ikC + ikwp C J2 + ip ( p C + Ew 2 k 4 /12 + wk 2 C /2 

13 * 2 55 32 *2 3 3 55 


(2.7) 


+ k 2 s 2 pw 3 /|2 + pws 2 ) 


where w = b/Z, p, E, are nondimensionalized quantities and p , p , , 



P 

3 2 


are defined in (1.11). 
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III. NUMERICAL INVERSION 

The inversions are accomplished by the Fast Fourier Transform (FFT) 
techniques [2] , which consists of a transformation from Laplace to Fourier 
transforms, and a two-dimensional numerical inversion using the usual 
FFT alogarithm. Notice the Laplace inversion formula 

C+i 00 

f(s) e St ds 

C-i°° 

r 

Set s = C + ia 

r 00 

f (C+ia) e Ct e iat da (3.1) 

-00 

r 

where C and a are both real, and C is greater than the largest real 
part of all singularities of f(s). Numerically, the double Fourier transform 
(or inversion) has the following form [1] : 


k 


f(t) = 


2ttx 


f(x,t). 


K K 

x t 

tt 2 NM 


- i[K x C 1 'l )x + K t (1 -^ t5 


N M 

Z Z f(I,J) e 
1=1 J=1 


(I-l)x (J-l) 

27TlL — n — + 


t] 


(3.2) 
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where N, M are the number of points in x and t direction respectively, 
and K , K are the corresponding half-frequency range. 

A t- 

For the present problem, the determination of C is through the following 
considerations : 

The form of inversion integrals are, in general. 


I 


1 

27Ti 


C+ ioo 


C-ioo 


FCk^s)- 
ACk i ,s) 



P = p(k^ ,s) , Re(p)>0 

(3.3) 


and A is given by (1.14). It is easy to see the singularities of the integral 
of (2.3) are: 

a) Poles at A = 0, 


A = (p -p ) ik (C -ty \p C ) + (i}j -i \) )(k 2 C +p p C ) = 0 

21 1 13 31 32 33 32 31 1 13 *1*2 33 


( 34 ) 


implies 


(P-P.)[C , k JP,P +c ) 2 + c ( c p 2 -s 2 -C k 2 )(C p 2 -s 2 -C k 2 )] 

21 1 3 1 1 2 1 3 5 5 3 3 5 5 1 1 1 1 5 5*2 11 1 J 


' ( C 1 , +C R )[p (C p-s-c k 2 ) - p (C p 2 -s 2 -C k 2 ) ] (k 2 C + p p C ) 

13 55 1 552 111 *2 55^1 111 113 ^1^233^ 


0 


( 3 . 5 ) 
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Interpretation of this condition is best understood 

for the particular case of an isotropic material, i.e. for 


It 


C = A + 2y = C , C = A, C = y, 

11 3313 55 


choose C = C = 

6 G 5 5 

is easy to see that from (1.10), (1.11) 


U 


det = [yp 2 - (yk 2 + ps 2 )] {(2y+A) p 2 - [(A+2y) k 2 + ps 2 ] } 

p 2 = (yk 2 + ps 2 )/y , p 2 = [ (A+2y) k 2 + ps 2 ]/(A+2y) 


* 


3 1 


= ik /p 

l r i 


ip = -p /ik , 

32 2 1 


1/2 1/2 
Here [y/p] = v s , is the shear wave speed and [(A+2y)/p] v^, is the 

longitudinal or pressure wave speed for isotropic materials. 
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Thus for the isotropic case that, A = 0 implies 


4pp 2 + 




P, 

+ — (X + 2y)] = 0 

k 2 


(3.6) 


where p^ p 2 are defined above. This is the equation for the Rayleigh wave 

speed v^ .= (is /k.) which is found as a real root of (3.6) (see e.g.[3]). 

For the case X = y (Poissons ratio = 0.25), v D = 0.919 v . 

R s 

For the anisotropic case a computational scheme to calculate the zeroes 
of (3.5) has been written, (Figure 2). 


b) Branch points : The branch points of the integrand in (3. 3) are the same 

as those of the functions p (k ,s), p (k ,s) , thev are: 

11 2 1 

i) p = 0, or p — 0 which implies (C k 2 + ps 2 ) (C k 2 + ps 2 ) = 0 
1 2 11 1 55 

i.e., k/s = ±i/p/C , ±i /p7'C pure imaginary. 

55 11 

These correspond to longitudinal and shear wave speeds for an isotropic material 


ii) p 


i >2 


= + fT± 


/ B 2 - 4AC /2A 


(3.7) 


C C > 0 

3 3 5 5 


with A = 
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b = - ps 2 (c +c ) + k 2 [ (c +c y - e c -c 2 ] 

3355 1355 113355 


C = (C k 2 + p s 2 ) (C k^ + Ps^) 

11 55 


2 . n „2- 


These branch points are those values of s/k which render B - 4AC = 0, 
and are branch points of second order. ^ The distribution of these points is 
shown in Figure 3 . It has been shown [5] that the contribution of these branch 
points to the value of the integrals (3.3) is important only when one considers 
the multi-reflected and refracted waves in layered media, or when the position 
of interest is very close to the impact origin. In this study we were more 
concerned with how the energy is propagated away from the impact point, which 
is mainly associated with surface waves, thus we ignored the contribution of 
these branch points.^ 

The contour of integration for the Laplace inversion is as shown in 
Figure 3. Notice the branch cuts are extended to negative infinity, in 
accordance with the requirement that C > max real part of the singularities. 
The requirement that real (p) 0 also determines the correct sheet 

of the Riemann surfaces. Numerically, since the branch points are located 

ct 

at s/k = constant, as k gets large C should be large, and the factor e 
in the Laplace inversion expression will rise sharply to an unmanageable size. 
Since, in the last paragraph, we have noticed the contribution of the branch 
points is unimportant, a path r 2 is chosen to replace Tj by the Cauchy's 
integral theorem. Notice the advantage of integration along f 2 is that C q 
is a positive constant independent of k. The determination of optimal 
is discussed in [2], Here, in order to minimize the aliasing and round off 
errors in numerical computations simultaneously, we choose 
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c = _J_ £n (g/r) 

° 3MX t 


where g = p /C , r = 10 x t— . A is & small ti 
6 V Q 66 A. t 


than the impact contact time. Further details are given in 


(3.9) 

interval, less 
Appendix A. 
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IV. RESULTS 

A computer program has been written to calculate the stresses in a 
plate with an elastic beam on one edge under a transient impact load 
distribution along the edge. A program description, flow charts, input 
data formats, and sample printout of the program are contained in the 
appendices to this report. In this section we will summarize some of 
the results obtained from this computer program. These results were 
calculated for an anisotropic plate with effective elastic constants 
of 55% graphite fiber/ epoxy matrix composite obtained from Reference 7 
and summarized in the Table. 

No- stri p case: 

The Rayleigh wave can be seen in the stress t^ on the edge as shown 
in Figures 4, 5 for a graphite fiber-epoxy composite for layup angles 0, -f 
15°. After the initial contact time, the stress is observed to propagate 
with little change at a speed near the calculated Rayleigh speed 
(Figure 2). This wave can be observed in the computed output in the 
space- time (x^, -t) plane Figure 7, as a band of non- zero values along 
a diagonal from the upper left to the lower right comer of the x^, -t 
plane. Caution is urged in using this program since spurious waves 
can enter the calculations due to the periodic nature of the finite 
numerical Fourier Transform. These spurious waves are data bands which 
lie along diagonals from upper right to lower left. In other words, 
only disturbances emanating from the impact source in the upper left 
corner of the -t plane of Figure 7 should be valid. A computer map 
of the space time history of the edge impact stress is shown in Figure 6. 

A contact time of 35 psec and contact length 2 cm was used in these cal- 


culations . 
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As is characteristic of surface wave effects , the stresses due to 
impulsive loading on the edge decrease with distance from the edge. This 
is shown in Figure 8 for two different layup angles. The normal stress 
t 3 2 appears to decrease to about 1/4 of its value on the surface at a 
depth equal to one half of the loading length a. The rate of decay from 
the edge depends on the layup angle. Another characteristic of edge 
impact is the development of tension in the normal stress under the impact 
point. This is shown in Figure 9. Thus while the compression part decreases 
with distance from the edge a tension tail d eve lopes in the wave. 


The effect of layup angle on the impact stresses can be seen in 
Figures 10, 11. The stress t^ 1 at the edge is larger than the impact 
pressure and decreases as the layup angle goes from 0° to +_ 45° (Figure 10). 

Below the surface or edge, the peak normal stress t^ at = 0 is 
a minimum for layup angles near +_ 30°, while the shear stress t^ 3 increases 
as the fiber angle goes from 0° to +_ 45° (Figure 10). 

An unexpected result is the shift of the maximum normal stress t^ to 
points off the impact axis x^ = 0 for layup angles greater than about 
30. A pronounced peak in t^ versus x^ beyond the impact pressure 
foot print, can be seen' in Figure 11 for +_ 45° layup angle. 

Impact protection strip case: 

The effect of bonding an edge impact protection strip to the half 
plane is shown in Figures 12-16. In Figures 12, 13, the increase in the 
thickness of a steel strip produces a decrease in the interface stresses 
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t i i , t 33 but creates an interface shear stress at the strip-composite 
interface. This shear stress reaches a maximum for strip thickness less 
than the impact footprint length and decreases for greater strip thick- 
nesses. Thus, if the strip is too thin, debonding can occur under im- 
pact due to induced interface shear . 

In Figure 14 one can see that increasing the strip thickness decreases 
the peak normal stress t^ and redistributes the load over a longer 
length under the strip. However, while the peak compression stress i s 
decreased by the strip, tension is created which could also produce 
debonding of the strip from the composite. 

For the no- strip case a Rayleigh wave was seen to propagate along 
the edge relatively unperturbed (Figures 4, 5). With the strip present, 
(Fig. 15) this wave becomes distorted as time increases. In fact, the 
beam- strip boundary conditions introduce dispersion in the edge waves 
which make the Rayleigh edge wave velocity dependent in the effective 
wave length of the disturbance. 

Finally in Figurelb shear stress distributions along the strip- 
composite interface are shown for a thickness near the shear peak 
(b/a = 0.25) and another for b/a = 2.0. In the latter case the shear 
is distributed over a larger length resulting in a lower peak stress. 

Also the strip delays the time of maximum shear from 1/2 to 3/4 T. 



Summary of Results 


An analytical-numerical method has been developed to solve the response 

of a composite plate with a bonded edge strip to in-plane impact type 

forces on the edge. Results of computer simulations reveal the following: 

1) Rayleigh edge waves can propagate away from the impact site with tension 
and compression up to values of the impact pressure, depending on layup angle 

2) Normal to the edge, the initial peak compression pulse decreases as 
it propagates into the plate but a tension tail develops as it 
propagates away from the impact site. 

3) The edge stress t^ under impact is decreased as the fiber layup angle 
goes from 0° to +_ 45°. 

4) Protection strips of thickness less than half the impact length can 
develop large interface shear under impact. 

5) The normal and edge stresses t^, t^ at the edge can be decreased 
significantly by protection strips of thickness greater than the half 
impact length. 
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TAELE I. - STRESS-STRAIN COEFFICIENTS FOR 55 PERCENT GRAPHITE 


FIBER-EPOXY MATRIX COMPOSITE 


[All constants to be multiplied by 10® psi; data obtained from ref. 7. ] 


0° Layup 

±15° Layup 


27.95 0.3957 0.3957 0 0 0 

1.170 0.4601 0 0 0 

1.170 0 0 0 

0.3552 0 0 

0.7197 0 

0.3552 



24.56 0.4000 1.986 000 

1.170 0.4558 000 

1.374 0 0 0 

0.3552 0 0 

2.310 0 

0.3552 



±30° Layup 

16.48 0.4118 5.167 000 

1.170 0.4400 0 0 0 

3.093 000 

0.3552 0 0 

5.491 0 



±45° Layup 

8.197 0.4279 6.758 0 0 0 

1.170 0.4279 000 

8.179 000 

0.3552 0 0 

7.082 0 



0.3552 


0.3552 
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APPENDIX A 

The Determination of Parameter C q 
C onsider the Laplace inversion of a function g(t) as 

g(s) e St ds (A.l) 

are real* We can change (A.l) into the 

Fourier inversion formula 

*» 

C t 

g (C o +iib) e iaJt diu = e ° x(t) (A. 2) 

which is then inverted by the Fast Fourier Transform technique. In the 
FFT scheme, a continuous function g(C Q +iuO is discretized and the 
infinite interval of integration is truncated. The error due to trunca- 
tion depends on each problem but doesn’t depend on C Q , thus in the 
determination of C q , we will assume the truncation error is negligible. 

It is shown [2] that the discretizing of the transform in one domain will 
cause aliasing error in. the other domain, e.g. sampling g at N points in a 

frequency interval, 0 < « < ft , will produce a transformed function x p (t) 
which is periodic and which differs substantially from x(t) for large 
enough t. For even x(t) this difference can be shown [2] to be given by, 

x p (t) * x(t) + x(t-T) 

for 0 < t < T/2 where T = N/fl. 


g(t) * 


C t 
o 


2tt 


stt) ■ ik 


rC +i°° 
1 o 


C - ico 
o 


and let s = C +ioJ, C and a) 
o ' o 
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It has been shown [2] also that the aliasing error is approximated 

.C (T-2t) 

by E (t) = e ° g(T-t) for the Laplace inversion (A. 2). 

a 

Notice the aliasing error is a decreasing function of C q . 

The other source of error is of course the round off error in compu- 

C o t 

tation. Since we multiply the resulting x(t) by e to get g(t), 
the rounding error is of the form 

C T 

Er(t) = e r(t) . 


The error bounds are then 


= [Max Ea[t) | = e 


-C (T-2t) 

Max ( g(T-t) 


C T 

£ = |Max Er(t) | = e ° Max|r(t) 

2 


0 < t < T 


Equating £ and e , the optimal C is then 
12 ° 


£n (Max g(T-t)/Max r(t)) 
T-t 


Chosing t = T/4, therefore 


C o = 3T £n(g/r) ’ 


g = Max g(T-t), r = Max r(t) 


Not 


N 


ice, empirically, r s 2 - 10 on single precision IBM 360 


systems . 
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APPENDIX B 


A FLOW CHART OF THE PROGRAM 


1. Main Program 
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APPENDIX C 

NOTES ON COMPUTER PROGRAM 

The choice of scales is very essential to the success of the present 
computational method. It is noticed that the accuracy depends on the 
number of points employed and the range of frequency spectra covered. 
Considering limitations of both computer storage and time, a time-space 
grid of 32 x 64 points was chosen for t > 0, x > 0. Thus, the non- 
dimensionalization of all equations and quantities are both necessary and 
important to the obtaining of meaningful data from the limited grid size. 
The numerical inaccuracies introduced have several origins: 

1. Theoretically, error has been introduced by the neglecting of the 

outer branch points contributions. It has been shown, for isotropic 

- 2 

cases, the contributions of these branch points behaved like r at 

large distance from a delta function loading at the origin r = 0. 

_1 

Compared with the r \ decreases of the contribution of the residue. 

Thus, for small r, or near the origin, the errors might be significant. 

An asymptotic form of the behavior of small r has been deduced for a 
simple delta loading at origin on an isotropic half-space [ ] . It is 

shown the error thus introduced is of the order of 5% maximum response 
in stress. 

2. The aliasing error introduced through the periodizing of the functions. 

3. The round-off error in Laplace inversion along with aliasing error, 
have been discussed in the determination- of C q (Appendix A). It is 
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found that error 3 is more serious of the two. Hence in calculations the 
maximum non-dimcnsionalized time should be restricted to below 6 or 8 
for reasonably good results. 

4. hrrors due to reflections at the boundary. Since the space grid is 
finite, it has been observed that whenever a wave hit the boundary of 
chosen space, a sizable numerical error will start propagating in, as if 
the wave were reflected from the boundary. The basic reason is due to 
the periodization of the space (x ) domain. In computation, this 
error should be avoided. To correctly determine the extent of the space 
(x^) domain, a priori, recognization of significant wave speed (at which 
most of the energy travels) is important. Usually an estimation will be 
sufficient. Then the proper nondimensionalizing constants can be 


chosen . 
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APPENDIX D 

PROGRAM INPUT AND OUTPUT 


A. Method: , 

FFT alogarithm for numerical inversion 


B. Input Data Cards: 

Card 1. NTEST, NSTRIP, NP (315) 

NTEST = 1 A test program for FFT (2-D) (beam) . 

0 The present program 
NSTRIP = 1 With strip 

0 No strip 
NP 


= 1 

Calculate: 

u 

l 

only 

2 


U , 

u 



1 

3 

3 


u , 

U , t 



1 

3 3 3 

4 


U , 

U , t , ,t 



1 

3 3 3 * ’ 11 

5 


V 

u , t , t 



3 3 3 1 1 

= 1,9, 

RHO, ANGLE 

(8E10,4) 


1 3 


CC(I=1,9) Material constants of composite, in the order 
C , C , C , C , C , C , C , C ,C . 

11 22 33 44 55 66 12 13 23 

RHO Density of composite (£/cm 3 ) 

ANGLE (degrees) lay-up angle 

Card 3. VEL , DM, El, ANU, DEN (8E10.4) 

VEL - Velocity of incoming particle m/sec 

DM - Diameter of impacting object cm 

El - Youngs Modulus of impacting object (psi) 

ANU - Poissons ratio of impacting object 
DEN - Density of impacting object (gm/cc) 


(psi) 
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Card 4. NSTRES , NK3, DX3 (215, F10,3) 

NSTRES = 1 

NK3 : Number of steps in direction 

DX3: Ax /£, step size in x direction 

3 3 

Photo elastic fringe order computer map in the x-t plane. 


RO, W, ES 

(3E10,4) 

RO: 

density of strip beam (g/cm 3 ) 

W: 

depth of strip beam (cm) 

ES:. 

Young's modulus of strip beam (psi) 


C. Output: 

1. Test problem: Appendis III 

2. Values of displacements 1 : 

2 : u 

3 

(stress) 3 : 

4 : t 

1 3 

5 : t 

i i 

3. Relative magnitude computer maps of displacements and stresses. 


r 

> f 


t/t* 
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FORTRAN IV G LEVEL 21 


MAIN DATE = 74186 11/28/2 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

G 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


******** PROGRAM TO CALCULATE STRESSES DUE TO EDGE IMPACT OF A PLATE* 

THIS PROGRAM CALCULATES THE ELASTIC RESPONSE OF AN 
ANISOTROPIC PLATE TO AN IN-PLANE EDGE IMPACT FORCE ON X3=0.0 
WHEN NSTRIP =1 , THE PROGRAM PLACES AN IMPACT PROTECTION STRIP 
,0R ELASTIC BEAM ON THE EDGE. THE IMPACT FORCE IS A HALF SINE 
FUNCTION IN TIME AND IS NON ZERO FOR 0<T<TO, (MICROSEC) . THE 
FORCE IS DISTRIBUTED ALONG THE EDGE AS P0(1-Xl**2), WHERE 
XI IS NORMALIZED BY THE HALF WIDTH OF THE IMPACT CONTACT 
LENGTH. 

THE METHOD EMPLOYS A FOURIER TRANSFORM IN THE EDGE DIRECTION 
XI, AND A LAPLACE TRANSFORM IN THE TIME DIMENSION. THE TRANSFORM 
OF THE FORCING FUNCTION IS GIVEN IN THE SUBROUTINE CALCIJL AND 
THE SOLUTION IS OBTAINED USING A 2-DIMENSIONAL FAST FOURIER 
INVERSION ROUTINE CALLED 'FOURT'. 

THE OUTPUT FOR A GIVEN DEPTH X 3 CONSISTS OF DISPLACEMENTS 
U 1 , U3, AND STRESSES T33,T11,T13, IN THE X1-TIME PLANE . 

THE STRESSES ARE NORMALIZED BY THE ELASTIC CONSTANT C(6,fi). 

IMPACT FORCE 
I 
I 
I 
I 

V I V 
V I V 
VI V 
V 


>>> * * * * 

******* 

♦ * * * 

«< STRIP 

c 


c 

> XI 

c 


C 

1 

c 

COMPOSITE 

C 

1 

c 


c 

1 

c 

HALF SPACE 

c 

1 

c 


c 

V 

c 


C 

X-3 


INPUT DATA 

NN ( 1) --2.0* MAX XI DISTANCE, NN(2)--MAX NO. OF TIMS UNITS 
NSTRIP=0 , NO STRIP, .. NSTRIP=1 , WITH STRIP 
NP=1,--,5, CALCULATES U 1 , U3 , T33 ,T1 1 , T 1 3 , IN THAT ORDER 
NP=6, CALCULATES DISPX. OF A BEAM ON AN ELASTIC FOUNDATION 
CC(9), ELASTIC CONSTANTS OP PLATE IN THE ORDER 

C11,C22,C33,C44,C55,C66,C23,C13,C12,1N PSJ 
RHO, DENSITY OF PLATE IN UNITS GM/CC 

ANGLE--LA TUP ANGLE OF COMPOSITE PLAT E, DEG FOR INFO ONLY 
VEL, VELOCITY OF INCOMING OBJECT METERS/SEC 
DM, DIAMETER OF IMPACTING OBJECT -CM. 

E 1 , A NU , YOUNG'S MODULUS AND POISSON'S RATIO FOP. IMPACTING BODY 

DEN, DENSITY OF IMPACTING BODY 

NSTRESS=1 

NK3, NO. OF DEPTHS X3,(FOR NK3=1,X3=0) 

DX3 , INCREMENT IN DEPTH X 3 (NORMA LIZED BY AO) 

RC, DENSITY OF PROTECTIVE STRIP GM/CC 


* This program has two extra cards to override the Hertz contact time and 

contact length calculation. Remove cards #51, 52. 



main 


DATE = 74186 


11/26/2 


FORTRAN IV G LEVEL 


21 


00C1 
00 02 
000 3 

0004 

0005 

0006 

0007 

0008 
0009 
00 10 
0011 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


W, THICKNESS TO WIDTH RATIO OF BEAM 
E5, YOUNG'S MODULUS FOR BEAM 

WL- -FOURIER WAVELENGTH (CM) OF THE ORDER OF AO 0? LESS 
WT--FOURIER WAVE PERIOD (SEC) OF THE ORDER OF TO OR LESS 
CHOICE OF WL, tf T DETERMINES DX, DT — DX=WL/2, DT=WT/2 
FG, TRANSFORM OF NORMALIZED FORCING FUNCTION F(X1)*G(T) 

THIS IS PROVIDED IN PROGRAM BUT CAN BE CHANGED BY THE USER 

OUTPUT DATA 

TC,TC CONTACT TIME (SEC,1.E-6 SEC ) (FROM HERTZ THEORY) 

AO, A, HALF THE IMPACT CONTACT LENGTH ,CM 

FO, MAX IMPACT FORCE PROM HERTZ THEORY , NEWTONS 

E, WAVE SPEED IN PLATE SQRT (C66/FHO) 

OR WAVE SPEED IND BEAM STRI P, SQRT ( El /D EN) UNITS CM/SEC 
CL=SQRT (Cl 1/RHO) LONGITUDINAL WAVE SPEED ALONG EDGE IN PLATE 
CS=SQRT (C55/RHO) , SHEAR SPEED ALONG EDGE OF PLATE, CM/SEC 
CR, RAYLEIGH WAVE SPEED ALONG FREE EDGE OF PLATE 

DX, DT SPACE TIME INCREMENTS IN Xl-T SPACE UNI TS--CM AND SEC 
DATA (I, J) , NORMALIZED TRANSFORM OF ONE OF DISPL. OR STRESSES 
--BEFORE CALL FOURT, AFTER CALL FOORT DATA IS A 2 DIM MATRIX OF 
DISPL. OR STRESSES IN Xl-T SPACE DEPENDING ON VALUE OF K 
IN THE LOOP 'DO 4 K=1,NP' 

U1, 03, NORMALIZED DISPL. IN PLANE OF PLATE (E . G. U1/A0) 
T33,T11,T13, NORM. STRESSES (E.G. T33/C66) 

-NOTE-- T33 ON X3 = Q SHOULD REPRODUCE THE FORCING FUNCTION 
WHEN THERE IS NO STRIP 

--NOTE-- AS A CHECK T13=0 ON X3=0 WHEN THERE IS NO STRIP 
THE FRINGE ORDER MAP PLOTS THE DTFF IN PRINCIPAL STRESSES AND 
MAY BE USED TO COMPARE WITH PHOTOELASTIC EXPERIMENTS OR TO 
LOOK FOR POINTS OF MAX IN PLANE SHEAR STRESSES 
******** ************ 

THIS PROGRAM HAS BEEN WRITTEN BY F.MOON AND C-K KANG UNDER 

A GRANT TO PRINCETON UNIVERSITY FROM THE NASA LEWIS RESEARCH 

LAB. 

******************** 

--NOTE TO THE USER— IN THE OUTPUT MAPS OF STRESSES OR DISPL., 
YOU WILL NOTICE EANDS OF SIMILAR NUMBERS RUNNING FROM TH E 
UPPER LEFT CORNER TO THE LOWER RIGHTCORNER — THESE APE WAVES 
WHICH EMINATE FROM THE IMPACT POINT — HOWEVER- WAVES RUNNING 
FROM RIGHT UPPER TO LEFT LOWER CORNER ARE SPURIOUS DUE TO THE 
DESCRFTENESS OF THE NUMERICAL FOURIER INVERSION PROGRAM 
--ALSO DATA FOR TIMES NEAR TMAX AT THE BOTTOM OF THE MAPS ARE 
USUALLY SPUPIOUS AND SHOULD NOT BE USED 

COMMON /MHC/ D, DKl 

COMMON /MC/DK2, FG1,T0,C1 1 ,C 1 3 ,C 33,C 55 , RHO , RO, W , ES , NSTRT P, AO 
DIMENSION DATA ( 1 28, 32) , MM (4 0) ,CC (9) 

DIMENSION NN (2) 

DIMENSION RDATA (40) 

DIMENSION CIDATA (40) 

DIMENSION FRNGE (64,32) 

DIMENSION T1 1 (64, 32) ,T 33 (64, 32) ,T 13 (64, 32) 

COMPLEX DATA , S 

COMPLEX Pi, P2, S1,S2,C1,C2,B,C,D,SI 
COMPLEX DK2,FG, SLAP, DKl 
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FORTRAN IV G LEVEL 21 SAIN DATE = 74186 11/28/2 


0012 


CCMPLEX*16 BS 

0013 


NN ( 1 ) = 1 28 

0014 


NN (2) =32 


c*** 

READ IN AND SPITE OUT DATA AND PARAMETERS 


C*** 

RHO , AO ,T0 MUST BE IN C.G.S. UNITS, PO MUST BE CONSISTENT WITH 

0015 


CALL INDUMP 

0016 


READ (5,102) NTEST,NSTRIP,NP 

0017 


.READ (5,100) CC,RHO, ANGLE 

0018 


WRITE (6,490) CC 

0019 


WRITE (6,491) RHO, ANGLE 

0020 


READ (5, 100) VEL,DH,E1 ,ANU,DEN 

002 1 


P0=CC(6) 

0022 


C11= CC(1) -CC<7)**2/CC{2) 

0023 


C33= CC (3) -CC (9) **2/CC (2) 

0024 


C 1 3= CC (8) -CC (7) *CC (9)/CC (2) 

0025 


C55= CC (5) 

0026 


PI= 3.14159265 

0027 


RHO= FHO/6. 895*1. E-4 

0028 


READ (5,101) NSTRES,NK3,DX3 

0029 


WRITE (6,504) NSTRES,NK3, DX3 

0030 


SI= CMPLX (0., 1.) 


c ******** CALCULATE THE IMPACT CONTACT TIME , RADIUS, AND PRESSURE 


c ******** BASED on hertz contact theory 

0031 


R=DM/2.0 

0032 


R=R* 1.E-2 

0033 


E1=E1*6895. 0 

0034 


E2=CC ( 2) *6895.0 

0035 


DEN = DEN* 1 . OE 3 

0036 


ABASS=4./3. *PI*R**3*DEN 

0037 


AK2 = 4./3. *SQRT (R) *E 1/ ( ( 1 . 0- A NU* ANU) *E1/E2) 

0038 


ALP=5./4. *A MASS *VEL*VEL/AK2 

0039 


ALF= (ALF) **0.4 

0040 


TC=2 . 943*A LF/VEL 

004 1 


F0=1. 14* AM ASS*VEL*VEL/ALF 

0042 


A=SQRT (ALF*R) 

0043 


A= 1 . 0E2* A 

0044 


TC=TC* 1 . 0E6 

0045 


WRITE (6, 710) VEL, TC, A,F0 


Q**** **************** ****************** 


c*** 

CALCULATE THE NON-DIMENSIONAL PARAMETERS 

0046 


BE= CC (6) 

0047 


E= SQRT (EE/RHO) 

0048 


RE= RHO 


c 

DEFINE TRANSFORM SPACE AND DISTANCE-TIME SPACE 


c*** 

UNIT DISTANCE —CM. 

0049 


AO= A 

00 50 


T0=TC* 1 . E-6 


c 

TEST CASE T0== 35E-6, A0=1 CM 

0051 


A 0=1.0 

0052 


T0=35. E-6 

0053 


AL=AC 


c *** 

UNIT TIME — SEC. 

0054 


TE= AI/E 


C*** 

SMALLEST WAVELENGTHS 

0055 


WL= AC/1 .5 

0056 


WT=TC/10. 0 
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FORTRAN 


0057 

0058 

0C59 
0060 
006 1 
0062 
0063 


0064 

0065 

0066 

0067 

0068 

0069 

0070 
007 1 

0072 

0073 

0074 

0075 

0076 

0077 
OC 78 

0079 

0080 
0081 
0082 
0033 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 


0099 

0100 
0101 
0102 

0103 

0104 


IV G t£VEL 21 SAIN DATE = 74186 11/28/2' 

C*** DIMENSIONS OF TRANSFORM SPACE 

XK=2.0*PI*AL/WL ; 

TK= 2. 0*PI *TE/WT 

C*** NCN DIMENS I ON ALI Z E CONSTANTS 
AQ= AO/AL 
T 0=T 0/TE 
DT=WT/2.0 
DX= WL/2 . 0 

WRITE (6,621) WT,WL,TK,XK 

C****** ** ****************************** 

C CALCULATE RAYLEIGH WAVE SPEED 

X1=0.0 
X2= C55 
N=0 

DO 90 1=1,100 
N=N* 1 

X= (XUX2) /2.0 
D 1 = (C55-X ) / (Cl 1 -X) 

D 1 = SQHT(D 1) 

D2= (Cl 3) **2-C33* (Cll-X) 

D3=C55*C33 
D3=SQRT (D 3) 

F=D 1 *D2/D 3 
F = X ♦ P 

I F ( F ) 81,82,83 

81 X1=X 

GO TO 92 

82 GO TO 91 

83 X2= X 

92 CONTINUE 
F1=F/C55 

DIF F=ABS (FI) -1 . OE-4 
IF(DIFF) 91,91,90 

90 CONTINUE 

91 CONTINUE 

WRITE (6,701) Cl 1,C33,C55,C13 

701 FORMAT (//, 1 0X ,6H Cll = ,E12.4,I>H C33 =,E12.4,6H C55 =,E12.4, 

16H C 1 3 = ,E 1 2. 4 ,/) 

CL=C1 1/RHO 
CL= SQRT (CL) 

CS=C55/RHO 
C5= SQRT (CS) 

CR= X/RHO 
CR=SCRT (CR) 

RDS=CR/CS 

WRITE (6,702) CL , CS,CR, RDS,N 

702 FORMAT (//, 10X, 1 2HLONG. SPEED =,E12.4, 13H SHEAR SPEED =,S12.4 
1 // , 1 OX , 1 6HRA YLEIGH SPEED =,E12.4,8K CR/CS = , P10 . 5 , 5X , I 5) 

C***************** ********************* 

C*** CONSTANTS FOR STRIP CASE 
READ (5,100) RO,W,ES 
WRIT E (6 , 521 ) RO, W , ES 
RO= RO/6. 895*1 . E-4 
IF (RO) 51,51,50 
50 CONTINUE 

E= SQRT (ES/RO) 
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0105 


WRITE (6,510) E 

C 106 

51 

CONTINUE 

0107 


RO= RO/RE 

0108 


ES= ES/EE 

0109 


W= W/AO/AL 


c** ****** ******** ********************** 

0110 

205 

CONTINUE 

01 1 1 


C11= Cl 1/EE 

0112 


Cl 3= Cl 3/EE 

0113 


C33- C33/EE 

0114 


C55= C55/EE 

0115 


P0= PO/EE 

0116 


RHO = RHO/RE 

0117 


FG1 = P0*T0 

0118 

206 

CONTINUE 

0119 


N- N N ( 1) 

0120 


n- N N ( 2 ) 


c*** 

CALCULATE THE LAPLACE INVERSION PARAMETER 

0121 


HT= WT/TE 

0122 


RH= PQ*WT*1 . E6 

0123 


CR= 2./3./M/WT*ALOG(RH) 

0124 


C LA P= CH 


c** * 

CALCULATE THE SECOND ORDER BRANCH POINTS 

0125 


A= RHO* (C33-C55) **2*RHO 

0126 


E= -2.* RHO* ((C3 3+C55) * (Cl 3*Cl 3+2 . *C1 3*C55 
1* (C 1 1 + C55 ) ) 

0127 


F= (Cl 3*Cl 3*2 . *C13*C55-C1 1 *C33) **2-4.*Cl1 

0128 


BS= 1. D0*E*E-1 .D0*4.*A*P 

0129 


D= CCSQRT(BS) 

0130 


P1= « 5/A* (-E*D) 

0131 


P2= F/P 1/A 

0132 


P1= CSQRT (PI) 

0133 


P2= CSQRT (P2) 

0134 


WRITE (6,509) 

C 1 35 


WRITE (6, 506) P1,P2 

Cl 36 

204 

WRITE (6,508) CLAP 

0137 


K 2= N/2+1 

0138 


M 2= M/2+1 

C 1 39 


IF (NTEST.EQ.1) GO TO 211 

0140 

211 

IF (NSTRIP.EQ. 1) WRITE (6,514) 

c*** 

GENERATE THE TRANSFORMED EXPRESSIONS 

0141 


REWIND 2 

0142 


DO 1 1= 1 , N 

0 143 


D K 1 - 2. *XK/N* (I-.5) -XK 

0144 


DO 1 J s 1,11 

0145 


DK2F=2.*TK/N*(J-.5> -TK 

0146 


SL AP= CMPLX (CLAP,DK2F) 

0147 


DK2" -SI*SLAP 

0148 


IF (NTEST.NE.1) GO TO 201 

0149 


FG= PI* (1 ,+C EXP (-SLAP) )/ ( SLAP* SLAP* PI*PI) 
1CSIN(DK1)/DK1* (1.+DK1*DK1/(PI*PI-DK1*DK1) 


c** * 

TEST FOR A STRING 

0150 


DATA (I,J) =-FG/(SLAP*SLAP+DK1*DK1) 

0151 


GO TO 1 

0152 

201 

CALL CALCUL (0) 

0153 

1 

CONTINUE 


1 1/28/2 


*C33*C55 
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0154 

0155 

0156 

0157 

0158 

0159 

0160 
0161 
0162 

0163 

0164 

0165 


0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 
0181 
0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 


IV G LEVEL 21 BAIN DATE - 74186 11/28/5 

IF (NTEST.NE.1) GO TO 202 
WRIT! (6, 507) 

CALL FOURT (DATA,NN,2, 1, 1,0) 

DO 15 KI- 1 , N2 
8= -SI*PI/N* (N-1) * <KI-1) 

DO 15 KJ- 1 , M 

C= -SI* PI/M* (B- 1) * (KJ- 1) 

T = PI/T K* (KJ-1 ) 

DATA (KI,KJ) = DAT A (KI,K J) *XK*TK/PI/PI/N/M*C EXP (B) *C.EXp(C) 

15 DAT A (KI , K J) - D AT A ( KI , KJ) *EX P (CLAP*T) 

WRITE (6, 506) ( (DATA (I, J) , J= 1,M), 1= 1,N2) 

GO TO 203 

£*#************************** 

C*** MAIN LOOP FOR DIFFERENT DEPTHS X3 
202 DO 6 K3- 1 , NK 3 
X 3= (K3- 1 ) *DX3 

C*** LOOP FOR CALCULATING DISPLACEMENTS AND STRESSES AT A GIVEN DEP'T’H 
DO 4 K= 1 ,NP 
WRITE (6,680) 

WRITE (6,501) X 3 
REWIND 2 
DO 7 1= 1 , N 
DO 7 J= 1 , H 

READ (2) Pi , P2,C1,C2,S1,S2,DK1,DK2 
B- CEXP ( -X3*P 1) 

C= CEXP (-X3*P2) 

C* * * DISPLACEMENTS 

IF (K.EQ. 1) DAT A (I , J) = Cl*B + C2*C 
IF (K.EQ.2) DAT A (I , J) - C 1*S 1 *B+C2*5 2*C 
C*** STRESSES 

IF (K.EQ.3) 

1DATA (I,J) = (Cl 3*SI*DKl-C33*P1*Sl) *C1*B+ (Cl 3*SI*DK1-C3 3*P2*S2) *C2*C 
IF (K.EQ.4) 

1 DAT A (I , J) - (C11*SI*DK1-C13*P1*S1) *C1*B+ (C 11* SI*QK 1-C 1 3*P2*S2) *C2*C 
IF (K.EQ.5) 

1DATA(I,J)= C55* ( (-P1+SI*DKl*Sl) *C1*B+ (- P2+ST*DK 1* S2) *C 2 *r) 

IF (K.LT.6) GO TO 7 

C*** DISPL. FOR A BEAM CN AN ELASTIC FOUNDATION 
EF=0. 1 

C*** FORCING FUNCTION 

FG= FI*FG1* (1 . + CEXP ( ~SI*DK2*T0) ) / (PI*PI -DK2*DK2*T0*T0 ) * 

14./(0K1*A0) **2* (CSIN (CK1*A0)/DK 1-A0*CCOS (DK 1* AO) ) 

D=W*CK1*DK1/12. 0* (ES*W*W*DK 1 *DK1 -RO*W* W*DK2*DK2) - RO* W*D K2*D K2 
1 +EF 

IF (K.EQ.6) 

1 DAT A (I, J) -FG/D 
7. CONTINUE 

CALL FOUR! (DATA ,NN , 2, 1, 1,0) 

DO 3 KI= 1 ,N 2 
B- -SI*PI/N* (N-1) * (KI-1) 

DO 3 K J= 1 , M 
C= -SI*PI/M* (M-1) * (KJ-1) 

T= PI/TK* (K J- 1 ) 

DATA (KI # KJ) - DATA (KI,KJ) *XK*TK/PI/PI/M/N*CEXP ( B) *CEXP (C) 

E= CL AP*T 

3 DAT A (KI, K J) = DATA (KI, KJ) *EX P (E) 
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c* ** 

CALCULATE THE PHOTOELASTIC FRINGE ORDER 

019 7 


DO 290 1=1,40 

0198 


DO 290 J= 1, 32 

0199 


IF (K, EQ. 3) T33 (I, J) =REAL (DATA (I, J) ) 

0 20 0 


IF (K.ECo 4) Til (I, J) =HEAL (DATA (I,J)) 

0201 


IF (K. EQ. 5) T13(I,J)=REAL(DATA(I,J)) 

0202 


IF (Ko £Q» 5) GO TO 280 

0203 


GO TO 290 

0204 

280 

CONTINUE 

0205 


FNGK= (T1 1 (I, J) -T33 (I, J) J **2*4.0*T13(I,J) *T1 

0206 


FBNGE (I, J) =SQ8T (FNGE) 

0207 

290 

CONTINUE 

0208 

214 

WRITE (6,511) K 


c*** 

FIND THE MAXIMUM VALUE 

0209 


RS = 1 . E-3 

0210 


DC 14 1*1,40 

021 1 


DO 14 J= 1,32 

C212 


S= DATA (I, J) 

0213 


TP= REAL (S) /RS 

0214 


IF {ABS (TP) .LT. 1 .) GO TO 14 

0215 


RS= REAL (S) 

0216 

14 

CONTINUE 

0217 


WRITE (6,516) RS 


C 

PRINT REAL PART OF DISPL. AND 

0218 


NIJ=N/4 

0219 


N IJ= 10 

0220 


DC 310 1*1, NIJ 

0221 


11=1-1 

0222 


WRITE (6,600) 11 

0223 


DO 300 J= 1 , 32 

0224 


SR=REAL (DATA (I, J) ) /RS 

0225 


SIHG*AIN AG (DATA (I, J) ) /RS 

0226 


RDATA(J) =SR 

0227 


CIDATA(J) = SIHG 

0228 

300 

CONTINUE 

0229 


WRITE (6,620) (PC ATA <L) ,L= 1 , H) 

0230 


BRITE (6,650) 

0231 


BRI TE (6, 620) (CtDATA(L) , L=1 , 1 1 ) 

0232 

310 

CONTINUE 

0233 


BRITE (6,680) 

0 2 34 


IF (K.EQ.1) WRITE (6, 630) 

0235 


IF (KoKQ.2) WRITE (6,631) 

0236 


IF (K„ EQ» 3) WRITE (6, 6 32) 

0237 


IF (Ko EQ • 4 ) WRITE (6 , 6 33) 

0238 


IF (KoEC.5) WRITE (6,634) 

0239 


IF (K.EQ.6) WRITE (6,6 35) 

0240 


WRITE (6 ,640) DX , DT 



PLOT THE RELATIVE VALUES 

0241 


DO 12 J= 1 ,M 

0242 


DO 13 1=1,40 

C243 


S= DATA (I, J) 

0244 

13 

MM (I)= REAL (S) /RS * 100 

0245 


WRITE (6,515) MM 

0246 

12 

CONTINUE 


c*** 

PLOT A MAP OF PHOTCELASTIC FRINGE ORDER 

0247 


IF ( Ko LT <» 5) GO TO 4 


STRESSES 


11/28/2 
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0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 
026 1 
0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 
0 280 
0281 
0282 

0283 

0284 

0285 

0286 

0287 

0288 
0289 


0290 

0291 

0292 

0293 


0294 

0295 

0296 

0297 


IV <3 LEVEL 21 


MAIN 


DATE = 74186 11/28/2 


RS= 1.E-3 
DO 312 1*1,20 
DO 312 J= 1 , 20 
SS=FRNGE (I, J) 

TP=SS/RS 

IF (ABS(TP) .LT. 1.) GO TO 312 
RS*SS 

312 CONTINUE 

WRITE (6 ,680) 

WRITE ( 6, 516) ES 
WRITE (6,690) 

DO 320 J* 1 , 32 
DO 315 1*1,40 
3S=FRNG E (I,J> 

315 mm (i) = ss/rs* ioo 
WRITE (6,515) mm 

320 CONTINUE 
4 CONTINUE 
6 CONTINUE 

100 FORMAT (8E10.4) 

101 FORMAT (2I5,3F10.3) 

102 FORMAT (315) 

490 FORMAT (5X , 55H ELASTIC CONSTANTS Cl 1 ,C22,C 33,C44,C55, C66 , C12, C13, C3 
1 2 ,/9E12 . 4 ,/) 

491 FORMAT (5X, 25H DENSITY OF PLATE GH/CC--, FI 0 . 4, 5X, 20H FIBER LAYUP AN 
1GLE — , Fl 0 . 4) 

500 PORHAT (2F 12. 4 , 4E1 5. 4) 

501 FORMAT (5H X3= F10.4) 

502 FORMAT (10H STRESSES I3/(8E15.7)) 

503 FORMAT (5X5H DATA, 15X10 HCONST ANTS , 20X 10HPARAHETERS /8E15.7) 

504 FORMAT (5X 8HN ST RESS* I3,5X4HNX3* I4,5X4HDX3* F10.4,/) 

505 FORMAT (14H SIMPLE POLES ) 

506 FORMAT (8E15.7) 

507 FORMAT ✓ (1 5H THIS IS A TEST ) 

5C8 FORMAT (30H LAPLACE INVERSION PARAMETER* F12.4) 

509 FORMAT (28 H SECOND ORDER BRANCH POINTS ) 

510 FORMAT (5X, 27H LONG. WAVE SPEED IN BEAM = ,F12.3,7H CM/SEC) 

511 FORMAT ( 14H DISPLACEMENTS 14) 

512 FORMAT (24H RAYLEIGH SPEED CS/CR+I ) 

513 FORMAT (8F15.8) 

514 FORMAT (/50X 10HWITH STRIP) 

515 FORMAT (2X, 4013) 

516 FORMAT (16H MAXIMUM VALUE* E15.7) 

520 FORMAT (/,5X,23H IMPACT PR ESS URE (PSI) = , El 2 . 4 , 5X , 25fl 1/2 CONTACT I.F 
1 NGT H (CM) = ,F12.4,5X,23H IMPACT DURATION (SEC) = ,E12.4,/) 

521 FORMAT {/, 5X, 18H DENSITY OF BEAM* , F 12 , 4 ,5X , 17H NORM. THICKNESS* , 
.1F12.4,5X,15H BEAM MODULUS* ,E12.4) 

600 FORMAT (/, 20X, 1 8HNORM ALIZED DIST.= ,14) 

620 FORMAT (1 1F10.5) 

621 FORMAT ( 5X, 20H WAVELEN (TIME-SEC) = ,E12 ,4,5X,20H WA VELE N (DI ST-C M )= 
1,F12.4,/,5X,23H MAX FREQ NO(NON DIM)* , Fl 2. 4, 5X, 2 3 h MAX WAVE NO <NC 
2N DIM)* , P12.4,/) 

630 FORMAT (30X,16H DISPLACEMENT U1,/) 

631 FORMAT (30X, 16H DISPLACEMENT U3,/) . 

632 FORMAT (30X,1 1H STRESS T33,/) 

633 FORMAT (30X,1 1H STRESS Til,/) 
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634 FOFMAT(30X,11H STRESS T13,/) /x 

635 FORMAT (30X,4lH TEST-DISPL . OF BEAM OR ELASTIC: FOONDATIOM,/) 

640 FORMAT(20X, 14H » Xl # DX= , F 1 2 . 4, / , 20X, 2H . l,/f 

1 2H |,/,20X,2H V,/,20X,2H V, /, 20X , 9HTI ME, DT- ,E12.4,//) 

650 FORM AT ( 1 0H I HAG PART) 

Ho FORM&T (/^20X # 17H FHINGE ORDER BAB 20X . 3 1H SQW ( (T1 1-T33) 

1 T 1 3*T 13)) 

710 FORMAT ( 5X,6H VEL= ,F12.4,5X,5H T0= ,F12.4,5X,4H A= ,F12.4,5X, 
15 H F0= , F 1 2. 4 ,/) 

203 STOP 
END 
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C0C1 

0002 

0003 

0004 

0005 

0006 

0007 

0008 
0009 
00 10 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 
0 04 8 

0049 
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C 


SUBROUTINE C ALCUL (NPOLE) 

COMMON /MMC/ D, DK1 

COMMON /MC/DK2,FGl,TO,Cl1,C13,C33,C55,RHO,RO r W,ES,NSTRtP,AO 
COMPLEX B1,B2 

COMPLEX GO,G01,HO,H01,G10,G11,H10 # H11 
COMPLEX 

COMPLEX DK1,DK2,P1,P2,S1,S2,C1,C2,SI,B,C,D , FG 
COMPLEX* 16 BS 
PI= 3.14159265 
SI= CMPLX (0 • , 1 . ) 

A- C33*C55 

B= C33* (RHO*DK2*DK2-C11*EK1*DK1) +C55* (RH0*DK2*DK2-C55*DK1*DK1) + 
1DK1*DK1* (Cl 3+C5 5) **2 

C- (RH0*D82*DK2-C1 1*DK1*DK1) * (RHO*DK2*DK2-C55*DK1 *DK 1 ) 

BS= 1. D0*B*B-1 .D0*4.*A*C 
D= CCSQBT(BS) 

Pi- • 5/ A* (- B + D) 

P2= C/PI/A 
P1= CSQRT(Pl) 

IF (PEAL(PI) .LT.0.O.AND.NPOLE.NE.2) Pl= -PI 
P2= CSQRT ( P 2) 

IF (REAL(P2) .LT. 0.0. AND. N POLE. NE. 2) P2=-P2 

S 1 - -SI/DK1 /Pi/ (C13+C55) * (RHO*DK2*DK2-C1 1 *DK1 *DR1 +C55 *P 1*P1) 

S2= -SI/DK1 /P2/ (C13+C55) *(RH0*DK2*DK2-C11*DK1*DK1+C55*P2*P2) 

IF (NPOLE.EQ.2) GO TO 201 
C*** FORCING FUNCTION 

FG= PI*FG 1 * (1. ♦CEXP(-SI*DK2*T0) ) / (PI*PI -DK2*DK2*T0*T0) * 

14./ (DK1*A0) **2* (CSIN (DKl*A0)/DK 1-AO*CCOS (DK 1*A0) ) 

201 IF (NSTRIP.EQ.1) GO TO 202 

D= (SI*DK 1*S 1-P 1) * (SI*DK1*C13-P2*C33*S2) - (ST*DK1*S2~P2) * (SI*DR1* 
1C13-P1*C33*S1) 

IF (NPOLE.EQ.2) RETURN 
206 Cl- FG/D* (SI*DK1*S2-P2) 

C2= FG/D* (Pl-SI*DK1*S1) 

GO TC 203 
202 CONTINUE 

C*** MATRIX FOR STRIP PROBLEM 

B1=H*W*(ES*DKl*DK1-RO*DK2*DK2) *DK1*DK1/12 . 0-RO*DK2*DK2 
E 1=W*B1 

B2=-N* (ES*DK1*DKl-RO*DK2*DK2) 

010- C55* (SI*DK1*S1-P1) 

H 1 0= C55* (SI*DK1 *S2-P2) 

G11=E2* (1. 0*SI*DK1*S1*N/2) 

H 1 1 = B2* (1.0+SI*DKl*S2*N/2) 

GO—- (SI*DK1*C1 3-PI *5l*C3 3) 

H0=- (SI*DK1*C13-P2*S2*C33) 

- G01=-SI*DK1*N*G 10/2.0 ♦ B1*S1 

H01--SI*DK1*H*Hl0/2. 0 ♦ B1*S2 
H 1 = B 10+H 1 1 
G 1 = G 10*G 1 1 
G=G0*GO1 


H=H0+H0 1 
D=G1C*H0-G0*H10 
D 

IF (NPOLE.EQ.2) 


♦ (H0*G11+H01*G10-G0*H1 T-G01*H10) *{H01*G11-H1 1*G01 
RETURN 


t 
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0050 

0051 

0052 

0053 

0054 


C1= -FG/D*H 1 
C2= FG/D*G1 

203 WRITE (2) Pi ,P2 ,Cl ,C2 ,S 1 , S2, DK 1, DK2 
RETUFN 

END . ■ 



LA 


ELASTIC CONSTANTS C 1 1 ,C22 , C3 3 ,C44 ,C55 ,C66 ,C 1 2 r C 1 3,C 32 , 

0.2456E C8 C.117CE 07 C.1374E 07 0. 3552E 06 0.2310E 07 0.3552E 06 0.4000F 06 'G.1986E 07 0. 4558E 06 

DENSITY OF PLATE GM/CC — 1.44CC FIBER LAYUP ANGLE — 15.0000 

NSTFESS= 1 NX 3= 1 DX3= . 1.0000 

VEL= 300.0000 T0= 87.3169* A= 1.3342* PO= 1042938.7500 

HA VELE N (TIME- SEC) = 0.3500E-05 HAVELEN (DIST-CM )= 0.6667 

MAX FBEQ NO (NON DIM)= 13.7654 MAX HAVE NO (NON DIH) = 9.4248 


C 1 1 = 0.2442E 08 C33 = C.1196E 07 C55 = 0.2310E 07 C13 = 0.1830E 07 


LCNG. SPEED = 0.1081E 07SHEAB SPEED = 0.3326E 06 

RAYLEIGH SPEED = 0.3020E 06 CR/CS = 0.90800 12 

DENSITY OF EE AM- 2.7000 NO RM . THICKNESS- 0.5000 BEAM MODULUS- 0.1124E 08 

LCNG. HAVE SPEED IN BEAM - 535757.750 CM/SEC 

SECOND ORDER BRANCH POINTS 

0.1997043E 02 0.0 0.0 0.2266004E 01 

LAPLACE INVERSION PAR AMETER= 0.5948 

HITH STRIP 

* This program has test cards (#51, 52) to override the Hertz calculation^ T0=3510~ 6 s A-l cm. 



STRESS T3 3 


-fc. 

Ov 


» XI, DX= 0.3333 

I 

I 

I 


V 

TIHE, DT= 0 . 1750E-05 


ooooooooooooooooo 

3 3 200000000000000 

14 12 720000000000000 

272315 6 0-1-1 0000000000 

41 36 24 11 2 -1 -2 - 1 0 0 0 0 0 0 0 0 0 

54 49 34 17 5 0 -3 -3 -2 -1 0 0 0 0 0 0 0 

67 61 44 23 8 0 -3 -4 -4 -3 -1 0 0 0 0 0 0 

78 71 52 29 12 2 -2 -5 -6 -5 -3 -2 0 0 1 1 0 

87 80 59 34 16 5 -1 -4 -6 -6 -5 -4 -2 0 0 1 1 

94 86 64 39 20 8 1 -3 -5 -7 -7 -6 -4 -3-1 0 1 

98 90 68 42 23 10 3 -1 -4 >6 -7 -7 -6 -5 -3 -1 0 

100 91 70 44 25 13 5 0 -3 -5 -6 -7 -7 -6 -5 -3 -2 

98 91 70 45 26 15 7 2 -1 -3 -5 -6 -7 -7 -7 -5 -4 

95 87 68 44 27 16 9 4 0 -1 -3 -5 -6 -7 -7 -7 -6 

89 82 64 43 27 17 10 6 2 0 -1 -3 -5 -6 -7-7 -7 

80 74 59 40 26 17 11 7 4 2 0 -1 -3 -5 -6 -7 -7 

70 65 52 36 25 17 12 8 6 3 1 0 -1 -3 -4 -6 -6 

57 53 43 31 22 16 12 9 7 5 3 1 0 -1 -3 -4 -5 

43 41 34 26 19 15 12 10 8 6 5 3 1 0 -1 -2 -4 

28 27 24 19 16 13 12 10 9 7 6 5 3 2 0 0 -2 

13 13 13 12 12 11 11 10 9 87 6 5 3 2 0 0 

1 13689 10 10 99876532 1 

-4 -3 -1 0 3689998876542 

-7 -6 -5 -3 0 2579998876 54 

-7 -7 -6 -4 -3 0 2 57 8 9 9 8 8 7 6 5 

-7 -7 -6 -6 -4 -2 0 2 4 7 8 9 9 8 8 7 6 

-6 -6 -6 -5 -5 -4 -2 0 1 4 6 8 9 9 9 8 7 

-6 -6 -6 -5 -5 -4 -4 -2 0 1 3 5 7 9 9 9 8 

-3 -3 -4 -4 -4 -4 -4 -3 -2 -1 0 3 5 7 9 9 9 

-5 -5 -4 -4 -4 -4 -4 -4 -3 -2 -1 0 2 4 6 8 9 

-1 -1 -2 -3 -3 -3 -4 -4 -4 -3 -3 -1 0 1 4 6 8 

-5 -4 -4 -3 -3 -3 -3 -3 -3 -3-3-3-2 0 1 3 5 


oooooooooooooooooo 
oooooo oooooooooooo 
oooooooooooooooooo 
ooooooooooooooo ooo 
oooooooooooooooooo 
oooooooooooooooooo 
oooooooooooooooooo 
oooooooooooooooooo 
ooooooooooooooooo o 
1 10000000000000000 
11 1100000000000000 
01 1210000000000000 
-2 0 0 1 2 11 00000000000 
-4 -3 -1 0 12210000000000 
-6 -5 -3 -1 0 1 2 2 2 1 0 0 0 0 0 0 0 0 

-7 -6 -5 -4 -2 0 1 2 2 2 1 0 0 0 0 0 0 0 

-7 -7 -6 -5 -4 -2 0 0 2 2 2 1 1 0 0 0 0 0 

-6 _ 7 _ 7 -7 -6 -4 -3 -1 0 1 22 2 1 0 0 0 0 

-5 -6 -7 -7 -7 -6 -5 -3 -2 0 1 2 2 2 1 0 0 0 

-3 -5 -6-7-7-7-6-5-4-200222 2 1 0 

-2 -3 -4 -6 -6 -7 -7 -7 -6 -4 -3 -1 0 1 2 2 2 1 

0 -i -3 -4 -5 -6 -7 -7 -7 -6 -5 -3 - 1 0 1 2 2 2 

1 0 -1 -2 -4 -5 -6 -7 -7 -7 -6 -5 -4 -2 0 0 2 2 

3 1 0 -1 -2 -3 -5 -6 -7 -7 -7 -7 -6 -4 -3-101 

4 3 2 0 0 -2 -3 -4 -6 -7 -7 -7 -7 -6 -5 -3 -1 0 

5 4 3 2 0 0 -1 -3 -4 -5 -6 -7 -7 -7 -6 -5 -4 -2 

6 5 4 3 2 1 0 -1 -2 -4 -5 -6 -7 -7 -7 -7 -6 -4 

7 6 5 5 4 3 1 0 0 -2 -3 -4 -6 -6 -7 -7 -7 -6 

8 7 6 6 5 4 3 2 0 0 -1 -3 -4 -5 -6 -7 -7 -7 

9 9 8 7 6 5 4 3 2 1 0 -1 -2 -4 -5 -6 -7 -7 

9 10 9 8 7 6 5 4 3 2 1 0 -1 -2 -3 -5 -6 -7 

7 9 10 9 9 7 6 5 4 4 3 2 0 0 -2 -3 -4 -5 



X3= 0.0 

DISPLACEMENTS 
MAXIMUM VALUE= 


3 

-0.7 1 04896E 00 


-p* 

“■4 


NORMALIZED DIST.= 0 


0.00050 

0.03979 0.14142 

0.27084 

0.41184 

0.54907 

0.67705 

0.78714 

0.67790 

0.94377 

0. 98577 

1.00000 

0.98948 0.95134 

0.89048 

C. 80461 

0.70045 

0.57593 

0.43953 

0.28894 

0.13502 

0.01142 

-0.04267 

-0.07263 -0.07276 

-0 .07738 

-0.06042 

-0.06357 

-0.03818 

-0.05201 

-0.01223 

-0.05377 


I BAG PART 

o.oooco 

-0.C000C 0.00000 

0.00000 

0.00000 

0.00001 

0.00001 

0.00002 

0.00001 

0.00002 

0.00003 

0.00052 

‘ NORMALIZED 
0.03504 0.12364 

DIST. = 
0.23925 

1 

0.36748 

0.49 373 

0.61244 

0. 71537 

0.80083 

0.86362 

0. 90441 

0.91968 

0.91198 0.87880 

0.82442 

0. 74692 

0.65225 

0.53870 

0.41383 

0.27579 

0.13423 

0.01894 

-0.03628 

-0.06768 -0.07048 

-0.07543 

-0.06060 

-0.06261 

-0.03975 

-0.05059 

-0.01542 

-0.04968 


XHAG PART 

0.00000 

-0.00000 O.COOOO 

0.00000 

0. 00000 

0.00001 

0.00001 

0.00002 

0.00001 

0.00002 

0.00003 

0.00046 

NORMALIZED 
0.02177 0.07726 

£IST.= 

0. 15621 

2 

0.24997 

0.34647 

0. 44004 

0. 52356 

0.59467 

0.64909 

0.68662 

0.70461 

0.70446 0.68451 

0.64749 

0.59235 

0.52308 

0.43886 

0.34489 

0.24037 

0. 13188 

0.03829 

-0.01809 

-0.05320 -0.06323 

-0.06962 

-0.06054 

-0.06012 

-0.04359 

-0.04736 

-0.02317 

-0.04079 

IflAG PART 

0.00000 

0.00000 0.00000 

O.OOCOC 

0.00000 

0.00000 

0.00001 

0.00001 

0.00001 

0.00002 

0.00002 

0.00028 

normalized 

0.00689 0.02688 

DIST . = 
0.06236 

3 

0.11326 

0.17235 

0.23427 

0.29336 

o'. 34651 

0. 39036 

0. 42373 

0.44479 

0.45370 0.44958 

0.43352 

0.40525 

0. 36662 

0. 31768 

0.26100 

0.19681 

0. 12822 

0. 06182 

0.00870 

-0.03042 -0.04999 

-0.06009 

-0.05857 

-0.05702 

-0.04745 

-0.04482 

-0.03114 

-0.03385 

I HAG PART 

o.ooocc 

O.OOOCC 0.00000 

O.OOOCO 

0.00000 

0.00000 

0.00000 

0.00001 

0. 0000 1 

0.00001 

0. 00001 

0.00030 

NORMALIZED 
0.00089 0.00200 

DIST • - 
0.00749 

4 

0.02457 

0.05280 

0.08844 

0. 12718 

0. 16551 

0.20047 

0. 23014 

0.25300 

0.26837 0.27564 

0.27493 

0.26619 

0.25006 

0.22685 

0.19759 

0.16280 

0. 12382 

0.08191 

0.03868 

-C. 00176 -0.03008 

-0.04672 

-0.05288 

-C. 05368 

-0. 04925 

-0.04457 

-0.03642 

-0.03159 

IMAG PART 

0.00000 

0.00000 -0. 000 GO 

-0.00000 

0. 00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00001 

0.00001 

0.00032 

normalized 

0 . CO 039 -0.00267 

EIST.= 

-0.01052 

5 

-0.01476 

-0.00909 

0. 00649 

0.02920 

0. 05590 

0.08350 

0. 10982 

0. 13298 

0.15198 0.16592 

0. 17462 

0.17771 

C. 17551 

0.16799 

0.15576 

0. 13895 

0. 1 1848 

0. 09467 

0.06546 

0.02911 -0.00428 

-0.02892 

-0. 04250 

-0.04904 

-0.04829 

-0.04595 

-0.03865 

-0.03469 

IMAG PAST 

0.00000 

0.00000 -0.00000 

-0. GOOCO 

-0.00000 

-0.00000 

-0.00000 

0.00000 

0.00000 

0.00000 

0.00009 
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X3- 0.0 

DISPLACEMENTS 4 

MAXIMUM VAIUE= -O'. 116906 IE 01 


0.00 258 

0.0059 1 

NORMALIZED 

-0.02645 

DIST. = 
-0.05231 

0 

-0.04372 

0.54124 

0.61882 

0.68126 

0.72641 


0.75286 

0.39153 

0.25986 

0. 15272 

0.07511 


0.01861 

I HAG PAST 

O.OOCCO 

0.00000 

0.00000 

-0.00000 


-0.00000 

0.00200 

0.00707 

NORMALIZED 

-0.00660 

DIST.* 

-0.01783 

1 

-0.00073 

0.59461 

0.66947 

0. 72778 

C. 76810 


0.78838 

0. 39149 

0.26578 

0.16154 

0.08 150 


0.02557 

IM AG PART 

C. 00000 

0.00000 

0. 00000 

-O.OOCOC 


-0.00000 

0.00095 

0.0 1745 

NORMALIZED 

0.05203 

E IST.= 
0.08323 

2 

0. 12706 

0.75508 

0.82171 

0.86766 

0.89339 


0.89515 

0.39022 

0. 28142 

0. 18630 

0. 10048 


0.04569 

IB AG PART 

0.00000 

0.00000 

0. 00000 

C.OOCOO 


0.00000 

0.00057 

0.03187 

NORMALIZED 
0. 10966 

DIST .= 
0.18732 

3 

0.26129 

0.89837 

0.95550 

0.98870 

1.00000 


0.98414 

0. 38626 

0.30196 

0.22148 

0.13213 


0.07587 

I HAG PART 

0 . ooocc 

0.00000 

O.OCOOO 

o.oocoo 


0.00000 

0. 0CC19 

NORMALIZED 
0.01996 0.08833 

DIST. = 
0.16952 

4 

0.24197 

0.76750 

0. 82278 

0.85892 

0.87669 


0.87127 

0. 38149 

0.32254 

0.25902 

0. 17532 


0.11222 

I HAG PART 

0.00000 

O.COCCC 

0. 000G0 

O.OCCCC 


C. 00000 

0.00069 

NCPMAilZFD 
0.00496 0.04046 

CIST. = 
0.09727 

5 

0. 15250 

0.51583 

0 . 5 7 1 C 9 

0.61642 

0.64875 


0.66607 

0. 38097 

0. 34064 

0. 29211 

0. 22 508 


0. 15284 

I MAG PART 

o.ooooc 

0.00000 

0.00000 

0.00000 


0.00000 


*- 

u> 


-0.00204 

0.75988 

-0.01623 

0.06720 

0.74715 

-0.04614 

0. 15474 
0.71530 
-0.05638 

0.25280 

0.66481 

-0.08178 

0. 35369 
0.59890 
-0.05544 

0. 45160 
0.51992 

-0.00000 

0.00000 

0. 00000 

0.00000 

0.00001 

0.00001 

0.04631 

0.78875 

-0.01399 

0. 11913 
0.76829 
-0. 04097 

0. 20892 
0.72861 
-0.05868 

0. 30838 
0.66949 
-0. 07437 

0.40951 

0.59492 

-0.06856 

0.50684 
0. 50860 

0.00000 

0.00000 

0. 00000 

0. 00000 

0.00001 

0.00002 

0.19085 

0.87552 

-0.00618 

0. 27497 
0.83186 
-0.02674 

0. 37173 
0.76861 
-0.06290 

0.47550 
0. 68362 
-0.05665 

0.57739 

0.58354 

-0.09664 

0.67297 
0. 48166 

0.00000 

0.00000 

0.00001 

0.00001 

0.00001 

0. 00002 

0.34028 
C.945Q3 
0.01 127 

0.43188 
0.88079 
-0 .00827 

0.53084 

0.79607 

-0.06039 

0.63438 

0.68776 

-0.04181 

0.73316 

0.56538 

-0.11182 

0.82397 

0.45448 

0.00000 

0.00001 

0.00001 

0.00001 

0.00002 

0.00003 

0.30671 
0. 84603 
0.04213 

0. 37709 
0.79712 
0.01049 

0.45380 

0.73005 

-0.04408 

0.53712 

0.64180 

-0.03756 

0.61997 

0.53914 

-0.10398 

0. 69932 
0.43862 

0.00000 

0.00001 

0.00001 

0.00001 

0.00002 

0.00002 

0. 19436 
0.66710 
0.08553 

0.23295 

0.65118 

0.03096 

0.27765 

0.61852 

-0.01399 

0.33095 
0. 56999 
-0.04177 

0.39123 
0 . 50604 
-0.07016 

0. 45436 
0.43337 

0.00000 

0.00001 

0.00001 

0.00001 

0.00001 
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JOB 0350 7/5/74 IBM 360-91 

EXECUTION TIME 34 SEC. 

25 PAGES 

CORE 150 K 




IMPACT PROJECTION '-IMPACT FORCE 

STRIP DISTRIBUTION 


FIGURE i 


GEOMETRY OF IN- PLANE EDGE IMPACT OF 
A COMPOSITE PLATE WITH PROJECTION STRIP 




FIGURE 2a 



LAYUP ANGLE 


FIGURE 2b 

Comparison of In-plane plate edge wave speed with body wave speeds 
fiber layup angle for SS% graphite fiber/epoxy matrix composite. 




- 1.0 


ON EDGE, X s -0 
0* FIBER LAYUP ANGLE 
NORMALIZED BY 5.42 Pd 
a* lcm.,To>35psec. 


FIBERS 


FIGURE 4 


Edge Stress Surface Wave for 0° Fiber layup Angle: 55 
graphite fiber/ epoxy matrix composite 





- 1.0 


-0.4 


t M vs X, 

ON EDGE, X^O 
±15° FIBER LAYUP ANGLE 
NORMALIZED BY 4.62 Po 
a« Icm., To *35 /i sec. 





V R = x3mm/^sec. 
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FIGURE 5 \ 

Edge Stress t n Surface Wave for ^ 15° Fiber Layup Angle: 5S% Graphite fiber/expoxy 


matrix composite. 



Computer Map of Edge Impact Pressure t Impact Time 87 psec 
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SPURIOUS DATA 


FIGURE 7 


Computer Map of Edge Stress t ^ Surface Wave in the Space-Time (x^, t) Domain 





FIGURE 9 

Stress t versus Time at Different Distances from the Edge 
under the Contact Point. 





FIGURE 10 


Matri™Composit e . VerSUS La ^P an * le ®S% Graphite Fiber/epox> 





X,/a 



NORMALIZED STRESSES IN COMPOSITE 
AT STRIP -EDGE INTERFACE 




DATA FOR 55% GRAPHITE FIBER /EPOXY MATRIX 
± 15° LAYUP - STEEL STRIP 

FIGURE 12a 

Effect of Edge Strip Thickness on Interface Stresses. 


NORMALIZED STRESSES IN COMP 

At strip -edge interface 



normalized beam thickness 


DATA FOR 55% GRAPHITE FIBER/EPOXY 'MATRIX 
± 45° LAYUP ANGLE - STEEL STRIP 

FIGURE 12b 

Effect of Edge Strip Thickness on Interface Stresses. 


1.0 


b/a 


1.5 


2.0 
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X 3 b 0 , a = I cm ,T 0 = 35/i.sec 


FIGURE 13 

Effect of Edge Strip Thickness on Interface Stresses. 
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t M /Po ON X s =0 

STEEL STRIP THICKNESS b/a=0. 
±15° FIBER LAYUP ANGLE 
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FIGURE 14 

Effect of Edge Strip Thickness on Stress 


^33 distribution Along the 





RAYLEIGH EDGE WAVE 

0° LAYUP ANGLE t„/Po, TIME t-1.3 To 



FIGURE 15 

Effect of Edge Strip Thickness on the Rayleigh Edge Wave Shape. 
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FIGURE 16 






Effect of Edge Strip Thickness of Interface 
Distribution Along the Edge. 

Shear Stress t ^ 








